function [varargout] = XuF(varargin)
%
% implementation eqs 1-4 of xu-f 2008///
% or you could say, of Dittman et al, 2000
% or Dittman and Regehr, 1998....
%
% call:
% XuF no arguments just runs with some default parameters
% XuF (mode, varargs):
% if mode = 1, the table is passed in the second argument, and the plot is
% generated.
% if mode = 2, the table is in the second, the train in the third and the
% recovery times in the 4th. For modeling/fitting
% if mode = 3, the table is in the second, the train in the third and we
% assume that ther eis no recovery data.
% P. Manis 6/2008.
%
colors = {'r','g','b', 'c', 'k'};
% fprintf(1, 'Xuf\n');
% define constants: (Table 1, Yang and Xu-Friedman, 2008, modified)
switch nargin
    case 0
        table.F = 0.4; % release probability (constant; no facilitiation)
        table.k0 = 1/1.75; % /s, baseline recovery rate from depletion (slow rate)
        table.kmax = 1/0.025; % /s, maximal recovery rate from depletion (fast rate)

        table.td = 0.05; %  time constant for calcium-dependent recovery
        table.kd =  0.7; % affinity of fast recovery proces for calcium sensor

        table.ts = 0.015; % decay time constant of glutatme clearance
        table.ks = 1000; % affinity of receptor desensitization for glutatmate
        % The large value means no desense occurs (otherwise, ks should be about
        % 0.6)

        table.kf = 0.6; % affinity of facilitiation process
        table.tf = 0.01; % make facilitation VERY slow

        table.dD = 0.02; % sets Ca that drives recovery(Ca influx per AP)
        % 0.02 yields rate-dep recovery in 100-300 Hz
        table.dF = 0.02; % sets Ca that drives facilitation

        table.glu = 0.3;

        freqlist = [100,200,300]; % in Hz
        modelist = 0; % if 0, use regular; otherwise use exponential.
        plotvsn = 1; % control plots

        traindur = 0.5; % sec
        rec = 0.001*[10,20,30,40,50,100,200,300,500,1000]; % recovery points
        newfigure('XuF', 'Xu-F');
    otherwise
        runmode = varargin{1};
        switch runmode
            case 1
                table = varargin{2};
                freqlist = [50,100,200,400]; % in Hz
                if(nargin > 2)
                    modelist = varargin{3};
                else
                    modelist = 0; % if 0, use regular; otherwise use exponential.
                end;
                plotvsn = 1; % control plots

                traindur = 0.5; % sec
                rec = 0.001*[10,20,30,40,50,100,200,300,500,1000]; % recovery points
                newfigure('XuF', 'Xu-F');

            case 2 % with recovery data

                table = varargin{2};
                xtb = varargin{3}; % cell array of times for stimulation.
                xtr = varargin{4}; % and recoveries (matched)
                freqlist = zeros(length(xtb), 1);
                for k = 1:length(xtb)
                    freqlist(k) = 1/mean(diff(xtb{k}));
                end;
                modelist = 2;
            case 3 % no recovery
                table = varargin{2};
                xtb = varargin{3}; % cell array of times for stimulation.
                xtr = {}; % and recoveries (matched)
                freqlist = (1:length(xtb))*100;
                modelist = 2;
            otherwise
                fprintf(1, 'XuF - mode not recognized\n');
                return;
        end;
end;

for mode = modelist
    n = 0;

    if(mode == 1)
        ntrial = 10;
    else
        ntrial = 1;
    end;

    for kfreq = 1:length(freqlist)
        freq = freqlist(kfreq);
        n = n + 1;
        switch mode
            case 0
                mark_face = 'w'; % regular are open symbols
            case 2
                mark_face = 'k';
            otherwise
                mark_face = colors{n}; % poisson are filled
        end;

        for nt = 1:ntrial % to allow averaging of poisson train trials
            switch(mode)
                case 0
                    pt = (1/freq)*ones(1, freq*traindur); % pulse train, msec.
                    pt = cumsum(pt);
                    pt = 0.010+pt-pt(1); % make starts all the same.
                case 1
                    pt=[];
                    pt(1) = randexpo(freq);
                    k = 2;
                    while (max(cumsum(pt)) < traindur)
                        pt(k)=randexpo(freq); %#ok<AGROW>
                        k = k + 1;
                    end;
                    pt=cumsum(pt); % that's the stimulus train.
                    pt = 0.010+pt-pt(1); % make starts all the same.
                case 2 % get from incoming data
                    pt = xtb{n};
                    if(~isempty(xtr))
                        rec = xtr{n};
                    else
                        rec = 0;
                    end;
                otherwise
            end


            esum = zeros(length(pt)+length(rec), 1);
            for j = 1:length(rec)
                ptt = [pt max(pt) + rec(j)]; % pick just one recovery pulse per train to prevent interactions

                D = ones(length(ptt), 1); % depletion (1-D)
                S = zeros(length(ptt), 1); % desens (1-S)
                E = zeros(length(ptt), 1); % EPSP
                glu = zeros(length(ptt), 1);
                CaDi = zeros(length(ptt), 1);
                CaFi = zeros(length(ptt), 1);
                Fn = zeros(length(ptt), 1);
                Fn(1) = table.F;
                glu(1) = table.glu;
                dt = 1; % just a very long time since the last pulse
                [D(1) CaDi(1), Fn(1), CaFi(1)] = dstep(dt, table, 1, 0.1,  table.F, 0.1);
                S(1) = table.ks/(table.ks + glu(1));
                E(1) = Fn(1) * D(1)*S(1);
                for i = 2:length(ptt)
                    dt = ptt(i)-ptt(i-1);
                    [D(i) CaDi(i), Fn(i), CaFi(i)] = dstep(dt, table, D(i-1), CaDi(i-1), Fn(i-1), CaFi(i-1));
                    glu(i) = (glu(i-1) + table.F * D(i) )* exp(-dt/table.ts);
                    S(i) = table.ks/(table.ks + glu(i));
                    E(i) = Fn(i) * D(i) * S(i);
                end;
                if(j == 1)
                    esum = E;
                else
                    esum(length(pt)+j) = E(end);
                end;
            end;
            if(nt == 1)
                ER = esum(length(pt)+1:end);
            else
                ER = ER + esum(length(pt)+1:end);

            end;
        end;
        ER = ER/ntrial;
        if(nargout == 0) % only plot if we are user commandline mode
            subplot(3,1,3);
            plot((max(pt) + rec), ER/esum(1), [colors{n} 's-'], 'markerfaceColor', mark_face, 'markersize', 1.5);
            set(gca, 'Ylim', [0 2.5]);
            hold on;
            plot(pt, E(1:end-1)/esum(1),  [colors{n} 's-'], 'markerfaceColor', mark_face, 'markersize', 1.5);
            ylabel(gca, 'E(normalized)');
            if(plotvsn)
                xp = 1:length(D);
            else
                xp = ptt;
            end;
            subplot(6,1,1);
            plot(xp, D, [colors{n} 'o-'], 'markerfaceColor', mark_face, 'markersize', 1.5);
            set(gca, 'Ylim', [0 1]);
            ylabel(gca, 'Frac Rel');
            hold on;
            subplot(6,1,2);
            plot(xp, Fn, [colors{n} 'o-'], 'markerfaceColor', mark_face, 'markersize', 1.5);
            set(gca, 'Ylim', [0 1]);
            ylabel(gca, 'Facil');
            hold on;
            subplot(6,1,3);
            plot(xp, CaFi, [colors{n} 'd-'], 'markerfaceColor', mark_face, 'markersize', 1.5);
            ylabel(gca, 'CaD (d) CaF (x)')
            %    set(gca, 'Ylim', [0 1]);
            hold on;
            plot(xp, CaDi, [colors{n} 'x-'], 'markerfaceColor', mark_face, 'markersize', 1.5);
            subplot(6,1,4);
            plot(xp, S, [colors{n} 'd-'], 'markerfaceColor', mark_face, 'markersize', 1.5);
            ylabel(gca, 'Desens(S)')
            %    set(gca, 'Ylim', [0 1]);
            hold on;
        end;
        xout{n} = [pt max(pt) + rec]; %#ok<AGROW>
        yout{n} = esum/esum(1); %#ok<AGROW>
    end;
end;
if(nargout > 0)
    varargout{1} = xout;
    varargout{2} = yout;
end;

% [CaDi'; CaFi']'
% this is the discreet step for solving the analytical equations as
% described in the papers. This version, including F, is from Dittman et
% al.
function [Dn, CaDn, Fn, CaFn] = dstep(dt, T, Di, CaDi, Fi, CaFi)
% calculate next D from Equations 15/16, Dittman et al. 2000
% dt is time since last step
% Di is D at last time step
% T is the table (parameters/constants).
% CaDi is the calcium term for depression
% CaFi is the calcium term for facilitiation

CaDi = CaDi + T.dD;
CaFi = CaFi + T.dF;
CaDn = CaDi*exp(-dt/T.td); % change with next step
CaFn = CaFi*exp(-dt/T.tf);

% r = 1/dt; % convert to hz (units are SECONDS).
eta = (T.kd/CaDi + 1)/(T.kd/CaDi + exp(-dt/T.td));
eta = eta^(-(T.kmax-T.k0)*T.td);
Dn = 1-(1-(1-Fi)*Di)*exp(-T.k0*dt)*eta;

Fn = T.F + (1-T.F)/(1+T.kf/CaFn);
